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Abstract 



We propose a new algorithm for solving a system of two nonlinear 

•^U ■ transcendental equations with two complex variables based on the Miiller 

^^ | algorithm. The two-dimensional Miiller algorithm is tested on systems of 

. . different type and is found to work comparably to Newton's method and 

Y? ' Broyden's method in many cases. The new algorithm is particularly useful 

in systems featuring the Heun functions whose complexity may make the 
already known algorithms not efficient enough or not working at all. In 
Cn ' those specific cases, the new algorithm gives distinctly better results than 

the other two methods. 

As an example for its application in physics, the new algorithm was 
used to find the quasi-normal modes (QNM) of Schwarzschild black hole 
described by the Regge- Wheeler equation. The numerical results obtained 
by our method are compared with the already published QNM frequencies 
l*i _ and are found to coincide to a great extent with them. Also discussed 

y-~ ■ are the QNM of the Kerr black hole, described by the Teukolsky Master 

equation. 



1 Overview 



jrt ' Solving a system of two complex- valued nonlinear transcendental equations nu- 

merically is a task with varying difficulty, depending on the non-linearity of 
the system, the types of functions involved and the dimension of the space 
determined by the system. There are many well-known iterative root-finding 
algorithms, but most of them are specialized and optimized to work with a 
narrow set of functions - for example polynomials or functions with real- valued 
roots. The two most heavily relied upon one-dimensional algorithms - the secant 
method and Newton's method (or the Newton-Raphson method, [TJ [3J 2] ) can 
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work with a wide set of complex valued functions under proper conditions (see 
0]), but they have their weak sides. Newton's method requires the evaluation 
of the function and its first derivative at each iteration. This increases the com- 
putational cost of the algorithm and it makes the algorithm unusable when the 
procedure evaluating the derivative of the function has numerical problems (for 
example see the discussion for the Heun functions below) or when derivative 
becomes zero or changes sign. The secant method avoids this limitation, but 
in the general case, it has lower order of convergence (~ 1.618) compared to 
that of Newton's method (= 2) and the convergence of both of them is strongly 
dependent on the initial guess. 

These problems of the algorithms are inherited by their multi-dimensional 
versions such as the generalized Newton- Raphson method ([3]) and the gener- 
alized secant method (Broyden's method, [5 ). Although those problems can 
have varying severity, there are systems in which those algorithms cannot be 
used effectively. There are also some novel approaches (see [6], [7]), but when 
they rely on the same one-dimensional algorithms, they are likely to share their 
weaknesses as well. Clearly there is a need for new algorithms that will enlarge 
the class of functions we are able to work with efficiently 

A great challenge in front of root-finding algorithms in modern physics can 
be found in systems including the Heun functions. The Heun functions are 
unique particular local solutions of a second-order linear ordinary differential 
equation from the Heun type [HI [SJ EH [H] which in the general case have 4 
regular singular points. Two or more of those regular singularities can coalesce 
into an irregular singularity leading to confluent differential equations and their 
solutions: confluent Heun function, biconflucnt Heun function, double confluent 
Heun function and triconfluent Heun function. The Heun functions general- 
ize the hypergeometric function (and also include the Lame function, Mathieu 
function and the spheroidal wave functions (TOlIll]) and their wide applications 
in physics ([11]) was summarized recently in [12 a . From that paper and the 
cited therein, it is clear that the Heun functions will be encountered more and 
more in modern physics from quantum mechanics to astrophysics, hence we 
need adequate numerical algorithms able to deal with them. 

The work with the Heun functions, however, is more than complicated. 
While there are analytical works on the Heun functions, they were largely ne- 
glected until recently and therefore the theory is far from complete. The only 
software package currently able to work with them is MAPLE. Although, the 
routines that evaluate them in that package work well in the general case, there 
are some peculiarities - there are values of the parameters where the routines 
break down leading to infinities or to numerical errors. The situation with the 
derivatives of the Heun functions is even worst - in some cases, they simply 
do not work or their precision is lower than that of the Heun function itself. 
Also, in some cases there are no convenient power-series representations and 
then the Heun functions are evaluated in maple using numerical integration. 
Therefore the procedure goes slowly in the complex domain (compared to the 
hypergeometric function) which means that the convergence of the root-finding 
algorithm is essential. 



Despite all the challenges in the numerical work with the Heun functions, 
they offer many opportunities to modern physics. For example, they occur in 
the problem of quasi-normal modes (QNM) of rotating and non-rotating black 
holes. In this case, one has to solve a two-dimensional connected spectral prob- 
lem with two complex equations in each of which one encounters the confluent 
Heun functions. The analytical theory of the confluent Heun function is more 
developed than that of the other types of Heun functions, but still many un- 
knowns remain. Again, the evaluation of the derivative of the confluent Heun 
function is problematic in MAPLE and the behavior of the function is hard to 
predict. In that situation Newton's method cannot be used as a root-finding 
algorithm. Broyden's algorithm works well in most cases, but it is slowly conver- 
gent even close to a root. It is clear, then, that we need a novel algorithm, that 
will offer quicker convergence than Broyden's algorithm, but without relying on 
derivatives. 

To solve this problem, in the case of a system of two equations in two vari- 
ables, our team developed a two-dimensional generalization of the Miiller algo- 
rithm. The one-dimensional Miiller algorithm ([13 ) is a quadratic generalization 
of the secant method, that works well in the case of a complex function of one 
variable. It has very good convergence for a large class of functions (~ 1.84) and 
it is very efficient when the starting point (the initial guess) is close to a root. It 
is also well convergent when working with special transcendental functions. The 
two-dimensional Miiller algorithm seems to inherit some of the advantages of its 
one-dimensional counterpart like good convergence and usability on large class 
of functions as our tests show. The new algorithm was used to solve the QNM 
problem in the case of a Schwarzschild black hole and it proved to work with- 
out significant deviations from the results published by Andersson ([H]) and 
Fiziev ([IS])- Also, preliminary results for the QNM of the Kerr black hole are 
discussed and for them we also obtain a very good coincidence with published 
results [IB] . 

The article is organized as follows: Section 2 reviews the one-dimensional 
Miiller algorithm and its two-dimensional generalization, in Section 3 we discuss 
some physical application of the method and the numerical results obtained 
with it and in Section 4 we summarize our results. In the Appendix, the new 
algorithm is tested on various additional and more simple examples to verify its 
functionality. 

2 The Miiller algorithm 

2.1 One-dimensional Miiller's algorithm 

The one-dimensional Miiller algorithm ([3] [13]) is iterative method which at 
each step evaluates the function at three points, builds the parabola crossing 
those points and finds the two points where that parabola crosses the x-axis. 
The next iteration is the the point farthest from the initial point. 

Explicitly, for every three points Xj-%, Xj-i, Xj and the corresponding values 



of the function f(x) — > fj—2, fj-i, fj, the next iteration Xj+i is: 
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is the root with bigger absolute value. 

The advantages of this algorithm are that it is simple for implementation, 
it works with complex values, it does not use derivatives and generally, it has 
higher convergence than the secant method. It also works very well with special 
functions such as the confluent Heun function. For example the spectra in [TS] 
and [17] was obtained by the authors using that method. 

We will indicate the one-dimensional Miiller algorithm by the map: 

It: n{x in , F{x)) -+ X fin , 

where x m , x^ m are the starting and end points of the algorithm and the integer 
P is the number of iterations in which the algorithm completes. 

In this notations, the exit-condition for the one-dimensional Miiller algorithm 
will be | Xj—Xj-i |< 10~ d , where d is the number of digits of precision we require. 
Combined with checks for the value of the function (if it is small enough), we 
found this to be the best exit-condition for a root-finding algorithm since it 
works independently of the actual numerical zero in use, which may vary for 
the confluent Heun function. 

2.2 Two-dimensional Miiller's algorithm 

The two-dimensional Miiller method comes as a natural extension of the one- 
dimensional Miiller method. 

For two complex- valued functions F\ (x, y) and F^{x, y) we want to find such 
pairs of complex numbers (x t , y 7 ) which are solutions of the system: 

(F 1 (x I ,y I ) = 
\F 2 (x I ,y I ) = 

where / = 1,... numbers the solution in use. From now on, we will omit 
the index /, considering that we work with one arbitrary particular solution. 
Finding all the solutions of a system is beyond the scope of this article. 

Consider the functions F\(x,y), F2(x,y) as two-dimensional complex sur- 
faces z — F\{x, y) and z — F%(x, y) in a three-dimensional space of the complex 
variables {x, y, z|^|. Normally, to solve the system, one expresses the relation 



1 Equivalently, we can consider four real surfaces in five-dimensional real hyperspace, which 
are defined by four real functions of four real variables 
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Figure 1: A block scheme of the two-dimensional Miiller algorithm. 
\(Ci,C<i,Ca) is the plane with equation z = C1X + C22/ + C3 that crosses trough 
the 3 pairs of points (xi, j/j) and the function F% evaluated in them. The plane 
v is defined by the equation 2 = 0. The one-dimensional Miiller algorithm, 
fi(t m , F(t)) — > tf %n , is applied on the function of one variable F(t) with starting 
point t m and final point t' m 



y{x) from one of the equations, then by substituting it in the other equation, 
one solves it for x and from y(x) one finds y. In the general case, however, this is 
not possible. The idea of our code is to approximately follow that procedure by 
finding an approximate linear relation y(x) between the two variables and then 
using it to find the root of function of one variable trough the one-dimensional 
Miiller algorithm. 

To find the linear relation y(x), at each iteration we form the plane passing 
trough three points of one of the functions and then the equation of the line of 
intersection between that plane and the plane z — is used as the approximate 
relation y(x). This basically means that the so found y(x) is an approximate 
solution of one of the equations which ideally should be near the real solution in 
the z = plane. Substituting this relation in the other function, we run the one- 
dimensional Miiller algorithm on it to fix the value of one of the variables, say 
x. Using the value of x in the first function, we again run the one-dimensional 
Miiller algorithm on it to fix the value of the other variable - y. Alternatively 
one can substitute the value of x directly in y(x) to obtain y. This ends one 
iteration of the algorithm. The process repeats until the difference between two 
consecutive iterations becomes smaller than certain pre-determined number. 
This process is systematized on Fig. (fTJ). 

Explicitly, the code starts by evaluating the two functions Fi^(xi,yi) in 



three starting pairs of points (i = 1, 2, 3) that ideally should be near one of the 
roots of the system. In our case, those three initial pairs are obtained from 
one starting pair to which we add and subtract certain small complex number. 
This artificial choice is done only in the first iteration (n = 3), afterwards we 
use the output of the last three iterations to form (x n —2,yn-2), (%n-ii Vn— i)j 
{x n ,yn) and the respective F\ y2 {x,y). Thus on every iteration after n = 3 the 
actual complex functions F\^ 2 {x n ,y n ) are evaluated only once outside of the 
one-dimensional Miiler subroutines. 

Next we construct the plane passing trough those three points for one of the 
functions, say F<i by solving the linear system: 

C\x n - 2 + C 2 y n -2 + C 3 = F 2 (x n - 2l y n - 2 ) 
C\x n -i + C 2 y n -i + C 3 = F 2 (x n -i,y n -i) 
C\x n + C 2 y n + C 3 = F 2 (x n ,y n ). 

From it one obtains the coefficients C\, C 2 , C3 of the plane z — C\x + C 2 y + C3. 

This plane is intersected with the plane z — (i.e. C\x + C 2 y + C3 = 0) and 
the equation of the line between those two planes is the approximate relation 
y(x) of the two variables. 

We substitute that relation in the first function Fi(x,y) — > Fi(x,y(x)) and 
we start the one-dimensional Muller on that "linearized" function of only one 
variable, x. After some pre-determined maximal number of iterations, the exit- 
ing point is chosen for x n+1 ( fi(x n ,Fi(x,y = y(x))) ->■ x n+ i 0) . 

Then, there are two possibilities. 

Algorithm Ml: one could use directly the relation y{x = x n +i) to find 
y = Un+i- Or, 

Algorithm M2: One can substitute x = x n +i in the other function F 2 (x, y) — > 
F 2 (x — x n+ i,y) in order to find y n +i using again the one-dimensional Miillcr 
algorithm (ji{y n , F 2 (x n+1 , y)) -^ y n+ i). 

Our numerical experiments showed that both approaches lead to convergent 
procedure. 

After (x n+ i, y n +i) are fixed, the two functions Fi y2 {x n+ i, y n +\) are evaluated 
and if the new points are not roots, the iterations continue. 

The exit-strategy in the two-dimensional Muller algorithm is as follows: 

1. To avoid hanging of the algorithm or its deviation to another root we fix 
maximal number of iterations for the one-dimensional Muller subroutine, 
P. Our experience shows that small P (3 -10) is often the more effective 
strategy, because it stops the one-dimensional Muller from straying too 
far from the actual root of the system. 

2. The precision-condition {\xj— Xj^i \< 10~ rf ) remains in force for the one- 
dimensional Muller. Usually the algorithm exits, because of j > P during 



2 Since the maximal number of iterations in the one-dimensional Muller algorithm is fixed, 
for simplicity we will omit the index P in this sub-section. The index of the iterations of the 
two-dimensional Muller algorithm is n. 



the first few iterations of the two-dimensional Miillcr and the closer to 
the roots it gets, the smaller number of iterations are needed in the one- 
dimensional Miillcr to reach d and to exit. 

3. The precision d defined by the absolute value of the difference between 
two consecutive pairs (x n ,y n ) (combined with the values of functions 
F 1 (x,y),F 2 (x,y) at them) - | x n - x n -i |< 10 rf , | y n - y„-i |< l(T d 
is the primary exit-condition of the two-dimensional Miillcr. When d be- 
comes smaller than certain value, the algorithm exits with a root. 

4. To avoid the two-dimensional Miillcr algorithm from hanging, we set a 
maximal number of iterations N after which the algorithm reaches exits 
without fixing a root. 

5. A common problem occurs when one of the functions becomes zero before 
the other function. In those cases, the algorithm accepts the fixed value 
for a root, say x^ m 1 and runs the one-dimensional Miillcr on the other 
variable until it fixes a root - /j,(y n ,F2(x^ m ,y)) — > yf m . The algorithm 
then exits with a possible root: (x^ m , yf m ). 

The procedure can be fine-tuned trough change in the starting pair of points, the 
initial deviation or by switching the places of the functions, or even by replacing 
the functions with their independent linear combinations. 

As wc will show in what follows, this method inherits some of the advantages 
of the one-dimensional Miillcr algorithm, like the quick convergence in proximity 
of the root and the vast class of functions that it can work with. The major dis- 
advantage comes from the complicated behavior of the two-dimensional complex 
surfaces defined by the functions Fi t 2{x,y) which require one to find the best 
combination of starting points and number of iterations in the one-dimensional 
Miillcr subroutine so that the algorithm converges to the required root (if it is 
known or suspected). Generally, it is hard to tell when one point is "close" to 
a root. In some cases, even if certain starting pair of points is close to a root in 
terms of some norm, using it as a starting point in the algorithm may still lead 
to convergence to another root or simply to require more iterations to reach the 
desired root than if other pair of starting points were used. 

It is important to note that unlike Broyden's algorithm and Newton's algo- 
rithm which are not dependent on the order of the equations in the system, our 
two-dimensional Miillcr algorithm depends on the order of the equations. The 
numerical experiments show that while for some systems, changing the places 
of the equations has little or no effect on the convergence, in other cases, it 
slows down or completely breaks down the convergence. While such inherent 
asymmetry certainly is a weakness of the algorithm, there are ways around it. 
For example, one may alternate the places of the equations at each iteration 
or use their independent linear combinations (F* 2 = ai^-Fi + /3i,2^2)- Those 
approaches make the algorithm more robust, but since they may cost speed, we 
prefer to set the order of the equations manually. 

A technical disadvantage is that the whole procedure is more CPU-expensive 
than Newton's method and Broyden's method, since it generally makes more 



evaluations of the functions - each one-dimensional Miiller makes at least 1 
iteration on every step of the two-dimensional Miiller, thus it makes at least 
4 evaluations of each function. This is because on each iteration of the two- 
dimensional Miiller algorithm the functions in use change and thus one cannot 
use previous evaluations to reduce time. Still, in some cases, as we will show, 
the so-constructed algorithm is quicker or comparable to Newton's or Broyden's 
method. 

3 Some applications of the method for systems 
featuring Heun functions 

We will work only with the confluent Heun functions, which are much better 
studied than the other types of Heun functions, due to their numerous physical 
applications. Besides their numerical implementation was used successfully in 
previous works by the authors. For details on the numerical testing, see the 
Appendix. 



3.1 First example 

As a first example, one considers the following system: 

F x {x,y) = UcunC{-1.3x,2y,l+x,4:X,l-ly-2x 2 ,.75y) = 

F 2 (x,y) = RcuTiC{9ix,2.3ix+y,2ix-l,-l.9x(i+y),2x 2 +2ix-l.3y-0.2,y) = 0, 

where HeunC is the confluent Heun function ( 8 ) in maple notations. 
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2.1991016319 + 0.2140611770J 
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14 
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80 
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92 
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Table 1: S numbers the system in use, t and N label the time and the iterations 

needed for the algorithms to exit. * denotes the roots dependent on the order 

of the equations in Ml and M2. 

The results for two pairs of starting points are presented on row S=l in 

table [TJ Here, the Newton method do not converge and cannot be used at all. 

This particularity is due to the behavior of the derivatives of the confluent Heun 

functions and the known problems of Newton's method near points where the 

derivative becomes zero or it is not continuous ([!]). From the table, one can see 

that in this case the Ml modification of the two-dimensional Miiller algorithm 



gives the best result, needing less than half the time of the Broyden algorithm 
to find a root. The M2 modification is slower, but still better than the Broyden 
algorithm ( t M1 < t M2 <t B ). 

We continue with two examples from the black hole physics: 

3.2 Physical Application — QNMs of non-rotating and ro- 
tating black holes 

The quasi-normal modes (QNMs) of a black hole (BH) are the complex fre- 
quencies (ui) that govern the late-time evolution of the perturbations of the BH 
metric ( PJ |2JB [2D [221 [23]). 

Second example: Rotating black holes 

To find the QNMs of a rotating black hole, one uses the exact solutions of the 
Teukolsky radial and angular equations, describing the linearized electromag- 
netic perturbations of the Kerr metric, in terms of confluent Hcun functions, as 
stated for the first time in full detail in [15] . From [24] , for the values of the 
parameters: s=-l, M=l/2, \r\ = 110, m=0, a=0.01, 9 = n/3, one obtains: 



Fi(sc,2/) =HeunC(-1.9996ix,2.0002ia;H-1.000Q0.0002ia;-1.000Q-1.9996a;(i+a;) ! 

1.9995x 2 + 1.9998ia; + 0.5000-y,-110.02e (4n24 " Qr9(:E)) + 1.0000) 
x (iio.oOe^- 7124 "^ 9 ^^ 2 - ^- 0002 '") =0 
HeunC' (0.04a;, -1.00, 1.00, -0.04a;, 0.50-1. 00j/+0.02x-0.0001x 2 , 0.25) 
2 ^ X,V ' ~ HeunC(0.04a;, -1.00, 1.00, -0.04a;, 0.50-1. 00y+0.02.x-0.0001a; 2 , 0.25) 
HeunC' (-0.04a;, 1.00, -1.00, 0.04a;, 0.50-1. 00y-0.02a;-0.0001a; 2 , 0.75) 



HeunC(-0.04a;, 1.00, -1.00, 0.04a;, 0.50-y - 0.02a;-0.0001a; 2 , 0.75) 



= 



For brevity, here the radial equation F\ (x, y) was rounded to only 4 digits 
of significance. In our numerical experiments, we used the complete system 
with software floating-point numbers set to 64, where the derivatives of the 
confluent Hcun functions HcunC were replaced with the associate Sn confluent 
Heun function according to equation (3.7) of [2"5"] . This was done to avoid the 
numerical evaluation HcunC so that the peculiarities of the numerical imple- 
mentation of the confluent Hcun function (i.e. the use of maple fdiff procedure) 
are minimized. The difference in the times needed to fix a root when HcunC 
is used and when it is not used is small for the modes (i.e. a;) with small imagi- 
nary part (At ~ 15s), but it increases with the mode number, until it becomes 
significant for modes with big imaginary part (for the 10 th mode - S=2.3 in 
table Q] - the difference is already At ~ 100s). This slowdown is due to the 
time-consuming numerical integration in the complex domain, needed for the 
evaluation of HeunC'. 

For that system, three pairs of starting points were used:(0.49 + 0.18i, 2.001 + 
O.li), (0.17 + 0.97*, 2.001 + O.li), (0.69e - 1 + 5.146i, 2.001 + 0.51e - li). The 
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results are in row S=2 in table Q] In this case, once again the two modifications 
of the two-dimensional Miiller algorithm Ml and M2 are much quicker than the 
Broyden algorithm (i M1 ~ t M2 < t B ). Newton's method do not converge to 
a root. The supremacy of the Miiller algorithms is clear and it is not isolated 
- it is observed for other modes or values of the parameters (for example, for 
m = 1). As for the precision of the method, the first two modes were found to 
coincide with at least 9 digits of significance with the already published results 
of electromagnetic QNMs of a Kerr black hole (see [16]). We could not find a 
published value for the third mode. 

Systems 1-2, Complex roots 



I * 

S=l | S=2 

I * 

* I ° 

I * 



Root number: 



Figure 2: A graphical comparison between the times needed by the Broyden 
method (denoted with asterisks) and the Miiller methods Ml (the solid circles) 
and M2 (the circles) 

Summarizing the results of the previous two examples (rows 1 and 2 in 
table [IJ in Fig. [5J one can see that the two-dimensional Miiller algorithm gives 
distinctly better times than the Broyden method. Newton's method in those 
cases cannot be used. 

Third example: Non-rotating black hole 

Finally, we consider in more detail the problem of the gravitational QNMs of 
a non-rotating black hole thus testing the new method on a very well studied 
physical problem. To find the QNMs, one uses the exact solutions of the Regge- 
Wheeler equations, describing the linearized perturbations of Schwarzschild 
metric, in terms of confluent Heun functions as stated for the first time in 
full detail in [T5] . From [T5] , when the mass of the BH is set to 2 ill = 1, we 
obtain the following system of the type (HJ: 



Fi = (cos(0)-l)(cos(0) + l)Le0eradreP(Z,2,cos(0))=O 
F 2 = HeunC(-2iuj, 2 iuj, 4,-2 u 2 , 4-Z-Z 2 + 2 w 2 , 1- |r | 



(2) 



-i((TT+e)/2±arg(u>)) 
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where a; is a complex frequency, I is the angular momentum of the perturbation, 
9 E [0,7r] is the angle which we set to 9 = it — 1CP 7 and \r\ — 20. 

The parameter e is a small variation (| e |< 1) of the phase condition arg(r)+ 
arg(w) = -^r/2 (defining the direction of steepest descent, see [15]) and we use it 
to study the stability of the so found QNM frequencies. Since the applicability 
of the second equation in the system Eqs.@, F 2 , may depend on the choice 
of the values of this parameter, the behavior of the solutions of Eqs.© under 
the variation of e is still an open problem. Below we present some preliminary 
results in this direction. 

Using Eqs. ((2|), we run the two-dimensional Muller algorithm to find the 
unknown I and uj. We set the number of digits that MAPLE uses when making 
calculations with software floating-point numbers to 32. The precision of the 
algorithms is still set to 15 digits. 

From the theory, it is known that I is an integer and I — 2,3..., thus, it offers 
a way to control the work of the algorithms. We looked for roots only for I = 2 
and we obtained 1.99(9) + 1 x 10~ 17 i, with the first different from 9 digit being 
the 17 th . 

Because the phase condition r = 



i(— 7r-fe)/2— iarg{u) 



r | e' v ^;/--~'»w includes the complex 
argument in non-analytical way, which cannot be differentiated, this problem 
cannot be solved directly using Newton's method. Broyden's method, however, 
also do not work well, because one of the roots y in the pair (x, y) is a real 
integer, while the other is complex. What happens is that the algorithm fixes 
the integer root very quickly and the finite differences which are used to evaluate 
the Jacobian become infinity. Because of this, the algorithm is able to fix only 
the first 10 — 11 digits, while the other algorithms fix 14 — 15 digits. Therefore, 
although Broyden's algorithm gives better times (see Table [2]) than the two- 
dimensional Muller algorithms due to the less evaluations of the confluent Heun 
function on each iteration, its precision is much lower and it cannot be increased 
even by raising the software floating-point number. 



Mode: 





1 


2 


3 


4 


5 


6 


7 


8 


9 


10 


t B [s] 


100 


99 


156 


196 


386 


240 


253 


282 


302 


368 


398 


<A/2 N 


317 


413 


595 


741 


1175 


799 


874 


892 


1364 


971 


1355 


tun N 


202 


218 


335 


357 


497 


457 


396 


613 


623 


594 


667 



Table 2: The times needed for Broyden's method (ts) and the two-dimensional 
Muller methods (£mi and £m2) to fix a root. Note that while the precision 
of the former is 10 digits, the precision of the other two is 14 — 15 digits. To 
obtain those times, we solve the system: [Fi + F 2 , F± — F 2 ] with starting points: 
U)\ri\ + 0.01 + O.Oli, 2.1 + O.Oli, where n = 0..10. 



The numerical results are summed in Table ((3]). In it, the QNM frequencies 
obtained from Sys. (J5J are compared to those found by Andersson ([33]) with the 
phase amplitude method. Recently, those results were confirmed by Fiziev (see 
[To] ) with the one-dimensional Muller method applied on the exact solutions of 
the radial equation in terms of the confluent Heun function for I — 2 . To check 
the accurateness of the new method, we evaluate A =1 u>Muiier2d — ui Andersson |- 
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n 


Our uj 


Andcrsson's w 


A 





0.7473433689+0. 17792463H* 


0.747343368+0. 177924630i 


1.68 x 10" 9 


1 


0.6934219938+0.547829750i* 


0.693421994+0. 547829714i 


3.60 x 10" 8 


2 


0.6021069092+0. 956553966i* 


0.602106910+0. 956553966i 


1.02 x 10" 9 


3 


0.5030099241+1. 410296405i* 


0.503009924+1. 410296404i 


1.01 x 10" 9 


4 


0.4150291596+1. 893689782i* 


0.415029160+1. 893689782i 


4.41 x 10~ 10 


5 


0.3385988064+2. 391216108i 


0.338598806+2. 391216108i 


9.67 x 10~ 10 


6 


0.2665046810+2. 895821253i 


0.266504680+2. 895821252i 


1.48 x 10" 9 


7 


0.1856446684+3. 407682345i 


0.185644672+3. 407682344i 


3.90 x 10" 9 


8 


-0.030649006+3. 996823690i 


0+3.998000i** 


0.0306 


9 


0.1265270180+4.605289542i 


0.126527010+4.605289530i 


1.44 x 10" 8 


10 


0.1531069502+5. 121653272i 


0.153106926+5. 121653234i 


4.52 x 10~ 8 



Table 3: A list of the frequencies we obtained for the QNMs of Schwarzschild 
black hole compared with the numbers found by Andersson. A =| LOMuiier2d — 
oJAndersson |- The first 5 frequencies (n = — 4, marked with *) were ob- 
tained also by Fiziev using the confluent Heun functions and coincide with the 
presented here except for the n — 1 where the published by Fiziev value is 
0.693421994+ 0.547829750i. The 8th mode, marked with **, was obtained by 
Leaver 28 . 



One can easily see that most of the modes coincide with precision more than 
10~ 8 . The mode with biggest deviation from the expected value is number 8 in 
the table [3] which is thought to correspond to the so called algebraically special 
mode. 

Algebraically special (AS) modes have a special place in the QNM studies. 
The Andersson method is not applicable for them and these are excluded from 
his consideration. Berti, Cardoso and Starinets ([MIE!]) make a review on the 
results so far concerning these modes. Theoretically the 8th mode should be 
purely imaginary with value 4z, if it indeed corresponds to the AS case. Leaver 
([28]) found that mode at frequency 0.000000 + 3.998000i. In our results the 
8th mode has small but non-zero real part. This clearly shows that mode does 
not agree with the theoretical predictions for the AS mode (see also [29]). 

Furthermore, the mode with n = 8 depends critically on the value of e 
in the phase condition. For e G [—0.1,0.1], our method could not find any 
corresponding frequency. Outside this interval, there is such mode, but the sign 
of its real part depends on the sign of e - positive epsilon leads to positive real 
part and vice versa (w„ =8 =sgn(e) 0.030649006 + 3.996823690i). 

We examine more closely the relation W„(e) for all modes: 

For modes n < 4 there is no dependence on e and they come in pairs sym- 
metrical to the imaginary axiqfp uj n = ±\Re(ui n )\ + Im(w n )i . These results 
confirm the roots for n = 0, 1, 2, 3, 4 published in 15 . 



3 The ranges where each mode is found, however, depend on e: for n = 0: e £ [+0.8, ±0.75], 
for n = 1 : e G [=F0.8, ±0.45], for n = 2 : e G [+0.8, ±0.25], for n = 3 : e G [+0.8, ±0.1], for 
n = 4 : e G [=F0.8,0], where the first sign corresponds to frequencies with positive real part, 
the second sign - to those with negative real parts. The imaginary parts for each mode n 
coincide. 
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Modes with n > 4 (but n / 8) depend on e and the symmetry is broken 
- there is only one frequency from each mode corresponding to certain e, i.e. 
w n( e ) = — sgn(e)|i?e(u; n )| + Im(ui n )i. The dependence of |i?e(w„)| and Im(u> n ) 
on e requires more careful investigation and is beyond the scope of this work. 
Outside the above mentioned ranges, the modes n > 4 disappear or get trans- 
lated. 

The so found relation w„ (e) may be a signal for some kind of instability that 
needs to be studied more carefully. The presented here results are preliminary 
and the detailed study will be published elsewhere. Here our goal is just to 
illustrate the ability of our new numerical algorithm to tread this problem. For 
the case n = 8, similar behavior was mentioned by other authors, too 26, 2"?ll2l)]. 

From the presented results, it is clear that the two-dimensional Miiller al- 
gorithms work very well for the complicated system of equations Eqs.([2]) and 
it even allowed the study of the dependency u;„(e) done here for the first time. 
The next step is to apply the method to the more complicated physical problem 
of QNM of rotating black hole (see Example 2.). As recently discussed by Fiziev 
(see [M|), rotating black holes are described by the radial and angular Teukolsky 
equations which can be solved analytically in terms of confluent Heun functions. 
Preliminary results show that our algorithm works well in this situation too. 

4 Conclusion 

We presented the general idea of a method for solving a system of two complex- 
valued nonlinear transcendental equations with complex roots based on the one- 
dimensional Miiller method. As it was shown, the new method inherits some 
of the advantages of its one-dimensional counterpart as good convergence for a 
large class of functions. The efficiency of the new method is illustrated on differ- 
ent examples (For examples including elementary functions and simple special 
functions see Appendix.). The numerical results showed that in many cases, it 
is comparable to Newton's method and Broyden's method, while in some spe- 
cific systems - like those featuring Heun functions - it is quicker than the other 
two. The complete mathematical investigation of the proposed new method, 
and especially its theoretical order of convergence under proper conditions on 
the class of functions F\ , F2 is still an open problem. 

Furthermore, the two-dimensional Miiller algorithm works well with special 
functions like the confluent Heun function in terms of which one can solve analyt- 
ically the Regge- Wheeler equation [To], the Zerilli equation [3D], the Teukolsky 
radial and angular equations [M] and many others. Using this algorithm, we 
were able to solve directly the problem of quasi- normal modes of a Schwarzschild 
black hole with higher precision than that of the Broyden method. The so found 
solutions agree to great extent with previous published numerical results thus 
confirming the usefulness of the method. The new algorithm was used to study 
the stability of the so found QNM with respect to small variations in the phase 
condition in the radial variable and some preliminary results were presented. 
For other applications of the method see the recent references [3T| [32] . 
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Future application of the two-dimensional Miiller method would be in the 
case of QNMs of a Kerr black hole described by the radial and angular Teukolsky 
equations, which will be published elsewhere. 
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Appendix 

5 Numerical testing 

All the algorithms are realized as procedures on the software package MAPLE, 
the tests are done on maple 15, on Linux x64, CPU Intel Centrino Core 2 Duo, 
on 2.2GHz. The number of digits that MAPLE uses when making calculations 
with software floating-point numbers is set to 64. For Newton's method and 
Broyden's method we used the analytical formulas [3], where the Jacobian in 
both cases is evaluated exactly or with finite differences respectively (i.e. with- 
out the Sherman-Morrison formula). 

The times presented below are obtained after running each procedure 10 
times using the garbage collection function gc(); in MAPLE on each calculation, so 
that each run represents an independent numerical experiment. The total time 
for each method is then divided by 10 and rounded to 3 digits of significance. 
This way, even though the times depend on the system load at the moment, they 
are representative for the four methods in each example. The notable exception 
of this "averaging" are all the systems featuring Heun functions, where such 
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procedure would require too much time and thus they are evaluated only once, 
using the function gc();. 

The precision in all the examples is 15 digits, but only the first 10 after the 
decimal point are presented here. In all the cases, the initial deviation where 
needed is 0.001. Some of the examples are from [18] p. 617-618. 

The numerical results for the test-systems are summarized in Table 2] in the 
Appendix. In it, we compare Newton's method, Broyden's method and the two 
versions of the two-dimensional Miiller algorithm discussed in section 2.2. 
M2 which uses two one-dimensional Miiller subroutines to fix the (x n +i, y n +i) 
and Ml which uses one one-dimensional Miiller subroutine to find x n +i and 
then it evaluates directly y n +i — (— C3 ~ Cix n +i)/C 2 . The number in the 
brackets in the Nmi and Nmi columns is P, the maximal iterations in the one- 
dimensional Miiller subroutine. Everywhere in the table, for each (x m , y m ), the 
four algorithms exit with the same (xf m ,y' ln ) with precision of 15 digits. 



5.1 Elementary functions 

1. F 1 {x,y) = y 2 +3x-5 + x 2 =0 
F 2 (x,y)=x 2 + 3y-l = 



2. F^x, y) = x(l - x) + Ay = 12 

F 2 (x,y) = (x-2) 2 + (2y-3) 2 = 25 



For those systems (rows S=l,2 in table HJ), the number of iterations and the 
time needed to find a root in the two-dimensional Miiller algorithms (for Ml in 
particular) are generally close to those of Broyden's algorithm (ttf <*B fa *Mi < 
tM2)- For real roots, however, the algorithms Ml and M2 are the quickest of 
the four. 

5.2 Trigonometric, exponential and logarithmic functions 

3. Fi(x,y) =y- l/4sin(x) - l/4cos(y) = 
F 2 (x,y) = 5x 2 -y 2 = 

4. F\ (x, y) — exp(—3x) cos(y) + x = 
F 2 (x, y) = x 2 - 3yx + y 2 = 



Again, for real roots, the two-dimensional Miiller methods are quicker than both 

other methods (rows S=3,4 in table EJ). Newton's method is much quicker when 

the initial conditions and the roots are complex. 

It is important to discuss the dependency of two-dimensional Miiller methods 
(Ml and M2) from the maximal number of iterations in the one-dimensional 
Miiller subroutine, P. In many cases, changing P only affects the time needed 
for the algorithms to complete. There are cases, however, where this parameter 
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becomes critical. For example, when one uses as starting points (4.4-5. Oi, 8.5- 
16i) on the system 5 = 5, the following 7 roots are obtained: 

r = (0.3487096094+0.4633971546i, 0.9129336096+ 1.2131895010J) 
n = (3.0248444374-4.3689275542i, 7.9191455477- 11.4380008313i) 
rf = (0. 1632674377+0. 6065137375i, 0.0623626119±0.2316676331i) 
r 3 = (1.1119158619- 1.8296636950J, 2.9110335191-4.7901217415J) 
r 4 = (4.0158133827- 5.6039287836J, 10.5135359284- 14.6712760260J) 
r 5 = (-5.3999170768-3.12 x 10~ 37 i, -14.1371664435-4.38 x 10~ 43 i) 

From them, Newton's method converges to the root r$ (after 25 iterations), 
Broyden's method - to r\ (after 27 iterations). With the two-dimensional Muller 
methods (Ml and M2) depending on the parameter P one obtains: 

• with M2: r+ for P = 3,9, 10, 12 - 14 and P > 16, r for P = 5, r 5 for 
P = 6, and r^ for P = 8 (after averagely 10 iterations), 

• with Ml: r 3 for P = 3, r 4 for P = 4-10 and P > 16 and r x for P = 11-14 
(after averagely 11 iterations). 

Similar behavior is observed in the next example: 



5. Fi{x, y) = ln(x 2 + y 2 ) - sin(yx) - ln(2) + ln(7r) = 
F2{x, y) — e x ~ v + cos(yx) = 0. 



For real starting points, Newton's method is not convergent, since it remains 
locked to the real axis. Broyden's method also did not exit with a root for real 
starting points. The two-dimensional Muller methods on the other side, when 
started from (0.5,0.5) gave the root rf (see below). Other purely real initial 
conditions cither gave a root or the algorithm did not converge. 

Using the complex initial conditions (2.27+0.001i, 1.27), the following four 
roots were found: 

r = (0. 2129109625-2. 4380400935i,-1.3216238026-4.6551486236i) 
rf= (0.9203224533±0.7487874838i, 1.4188731053+0.5453380689*) 
r 2 = (-1.6645201248+1.380553001i, 1.66452012482+ 1.38055300197i) 

Newton's method exits with Tq after 24 iterations. Broyden's algorithm is 
not convergent. From the two-dimensional Muller algorithms one obtains: 

• with M2 after averagely 12 iterations: rj~ for P — 3, r% for P = 4 and r^ 
for P > 5 . 



• 



with Ml after averagely 14 iterations: r^ for P = 3 and P > 5 and r x for 
P = 4. 



From the last two systems it is clear that P represents an additional pa- 
rameter of the two-dimensional Muller algorithm. It can be used to improve 
convergence, but in some cases, it can lead to different roots for the same initial 
conditions. Such instability depends on the system and it can be avoided by 
starting the algorithm closer to the root. 
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5.3 Special functions 

Finally, we consider the following two systems: 

6. F\ (x, y) = x 2 - y + 5 sin(a; - 2) = 
F 2 (x,y)=3(3,y) + 5x-3 = 



7. Fi(x,j/) =x 7 -e v - 
F 2 (x,y)=H 1 (7,y 



2 F 2 ([l},[3],x 2 -3x) = 
1 - x) = 



where J() is the Bessel functions of the first kind, H 1 ^) is the Hankel function 
of the first kind and 2 F 2 is the generalized hypergeometric function. 

In this case (rows S=6,7 in table H} the two-dimensional Miiller algorithms 
are comparable to Newton's algorithm, while Broyden's algorithm is often the 
quickest of the four. This is likely due to the computational burden of the deriva- 
tive or of each additional function evaluation. Note, however, that our goal is 
not to have an algorithm that is better than Newton's method, but to have an 
algorithm that has good convergence and that does not need to evaluate deriva- 
tives. In that, the performance of the new algorithms is satisfactory, especially 
since in some cases like 6.2 and 6.3, Miiller 's algorithms are the quickest. 

5.4 Discussion 
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Systems 1-7, Complex roots 
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Figure 3: A graphical comparison between the times needed by the Newton 
method, the Broyden method and the Miiller methods Ml and M2. With crosses 
we denote Newton's method, with asterisks - Broyden's method, with solid 
circles - Ml, with circles - M2 

The numerical investigations above (see Fig|3] and also table 0]) show that 
in general, the two modifications of the two-dimensional Miller algorithm work 
comparably well to the more established algorithms - Newton's and Broyden's, 
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even if sometimes they require more time and iterations to fix a root. In some 
specific cases, like those with real roots or those featuring confluent Heun func- 
tions (see Section 3), however, the two-dimensional Miiler algorithms are often 
the quickest of the four. 

It was demonstrated that while in Newton's and Broyden's algorithms the 
exit points depend only on the starting points (when the initial deviation is 
fixed), the two-dimensional Miiller algorithm depends also on the number of 
iterations in the one-dimensional Miiller subroutine (P). Surprisingly, the time 
for fixing a root do not depend in a straight-forward way from P, since sometimes 
increasing P leads to decreasing of the total time. This can be expected, because 
P is the maximal number of iterations in the one-dimensional Miiller subroutine, 
but the actual number of iterations depends on the precision. Therefore, this is 
one more way to fine-tune the algorithm to achieve a known or suspected root. 

The examples also showed the importance of the order of the functions in 
the system (see the systems in the table |4] marked with * and f). While in most 
cases all the four algorithms find at least one root of the system, there are initial 
conditions which lead to a divergence or to "undesired" root. This problem can 
be avoided by starting the procedure closer to the root, changing the places of 
the equations in the system or using a linear combination of the functions (say 
Fi—tFi+F2, F2— >F\— F^)- We conclude that even though the new algorithms 
admit some further improvements and developments, they work well enough to 
be tested in real physics problem. 
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s 


(„ \initial 


( x \ final 
/ {final 




t D j [si 

Broyden L J 

N 

"Broken 


*M2 [s] 
Nm2 


<MlN 

N M i 


1 


1.689 
-0.637 


1.1890465736 
-0.1379439181 


0.019 

8 


0.026 
10 


0.018 
10*(3) 


0.013 

8(3) 


1.321 + 3. 520i 
3.738- 1.927J 


0. 8214691720+3. 5201983985i 
4.2389950548 - 1 1 .9278229759* 


0.021 

8 


0.033 
10 


0.066 

12(3) 


0.033 

8(3) 


1.321-3. 520i 
3.738 + 1. 927i 


0. 8214691720-3. 5201983985J 
4.2389950548 + 1. 9278229759i 


0.021 
8 


0.033 
10 


0.063 

12(3) 


0.033 

8(3) 


2 


-.5 
3 


-1.0000000000 
3.5000000000 


0.020 

8 


0.032 
10 


0.019 

9(4) 


0.013 

9(3) 


3.046 
3.484 


2.5469464699 
3.9849974627 


0.019 

8 


0.035 
10 


0.020 
10(3) 


0.016 
9(3) 


0.726+4.335J 
-2.242-0.592i 


0.2265267650+4.3352949767J 
-1. 7424987313-0. 5927935709J 


0.024 
8 


0.037 
9 


0.033 

7(6) 


0.036 

8(6) 


3 


0.621 
-0.228 


0.1212419114 
0.2711051557 


0.039 
10f 


0.044 
13f 


0.023 

9(3)f 


0.019 

8(4)*f 


-0.422 + 1.476J 
-2.562+3.30H 


-.9222203725 + 1 .4764038337J 
-2.062147443+3.3013393343i 


0.045 
9 


0.061 
11 


0.081 
8(4) 


0.086 

11*(3) 


1.468- 1.635J 
-2.665+3. 656i 


0.9685241736- 1.6351708695i 
-2. 1656858901 +3. 6563532190J 


0.043 
9 


0.062 
11 


0.077 
7(5) 


0.086 

11*(3) 


4 


-.35 
-1.05 


-0.5600551872 
-1.4662435158 


0.045 
11 


0.047 

12 


0.032 

7(4) 


0.035 

10(4) 


0.55 - 0.6i 
1.14 -H 


0. 3487096094-0. 4633971546J 
0.9129336096-1.213189501i 


0.048 
10 


0.067 

12 


0.082 

7(6) 


0.097 

13* (3) 


6 


1.2 + 0.09J 
-5.5 + 0.01? 


.6863031247 
-4.3646459533 


0.141 
9f 


0.134 
Ht 


0.220 
9*(4)f 


0.271 

ll(3)f 


7.2 - 3.6i 
-11.9 + 5.00H 


5. 8404591703-3. 0854927956J 
-10.6712592035 + 5. 7445552813J 


0.188 
11 


0.171 

14 


0.156 
10*(3) 


0.198 

14* (4) 


-5.1 - 1.006J 
16.0 + 5.51i 


-4.9297777922 - 1 . 1922443124i 
17.4620338366 + 5. 7870418188J 


0.196 
11 


0.191 

15 


0.161 
11(3) 


0.168 
13(3) 


7 


1.1 - .45i 
-2.4-4.2i 


. 8288091 244 - .4046494664J 
-2.3507488745-4.6830120304J 


0.476 
10 


0.249 
13 


0.640 

11(3) 


0.333 

12(3) 


.5 - .87i 
-3.21 -5.14i 


0. 2656154750-0. 8757700972J 
-2. 9139425238-5. 1541326612i 


0.457 
10 


0.233 
13 


0.577 

8* (4) 


0.316 

11*(3) 



Table 4: S numbers the system in use, t and N label the time and the iterations 
needed for the algorithms to exit. * denotes the roots dependent on the order 
of the equations in Ml and M2. In the f case, the places of the equations were 
switched to obtain that root. 
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